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Active gels are a class of biologically-relevant material containing embedded agents that spontaneously generate forces acting 
■on a sparse filament network. In vitro experiments of protein filaments and molecular motors have revealed a range of non- 
equilibrium pattern formation resulting from motor motion along filament tracks, and there are a number of hydrodynamic models 
purporting to describe such systems. Here we present results of extensive simulations designed to elucidate the microscopic basis 
underpinning macroscopic flow in active gels. Our numerical scheme includes thermal fluctuations in filament positions, excluded 
volume interactions, and filament elasticity in the form of bending and stretching modes. Motors are represented individually 
as bipolar springs governed by rate-based rules for attachment, detachment and unidirectional motion of motor heads along the 
filament contour. We systematically vary motor density and speed, and uncover parameter regions corresponding to unusual 
statics and dynamics which overlap but do not coincide. The anomalous statics arise at high motor densities and take the form 
of end-bound localized filament bundles for rapid motors, and extended clusters exhibiting enhanced small-wavenumber density 
fluctuations and power-law cluster-size distributions for slow, processive motors. Anomalous dynamics arise for slow, processive 
motors over a range of motor densities, and are most evident as superdiffusive mass transport, which we argue is the consequence 
of a form of effective self -propulsion resulting from the polar coupling between motors and filaments. 



.1 Introduction 

Living matter fundamentally differs from dead (or passive) 
media in that it is driven by spontaneously-activating inter- 
nal processes, depleting some energy reservoir maintained 
by the organism's metabolism- ~—. This should be contrasted 
with passive materials which may be driven externally by e.g. 
an imposed boundary, or simply agitated by thermal noise. 
An immediate and crucial consequence with regards quanti- 
tative modeling of living systems is that they do not obey the 
principles of thermodynamic equilibrium^, necessitating the 
development of novel analytical and theoretical principles if 
such systems are to be understood on the same level as equi- 
librium matter. Ideally, such a program should proceed via in- 
cremental improvements to theory in tandem with experimen- 
tal verification and guidance, but the enormous complexity of 
real organisms eliminates facile comparison with any suitably 
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transparent theory. It is therefore often prudent to treat in vitro 
systems of known composition and reduced complexity. 

This approach has been applied with some success to the 
cellular cytoskeleton, a dynamic scaffolding of protein fila- 
ments and associated proteins that contributes to the mechan- 
ical, structural and motility properties of eukaryotic cells-^"— . 
This can be classified as an active gel, both because it contains 
molecular motors that generate (pN) forces on the filaments, 
and because the filaments themselves can translate or 'tread- 
mill' due to different growth rates at either end. It is possible 
with reconstituted in vitro active gels to inhibit treadmilling 
and add permanent biotin-avidin crosslinks between filaments, 
resulting in a static gel with both thermal and athermal sources 
of noise, the latter deriving from the action of motor proteins 
on the network^^"— . These athermal noise sources have been 
shown to violate the fluctuation-dissipation relation^ii, cate- 
gorically placing active gels outside the realm of equilibrium 
thermodynamics, and the athermal noise spectrum can be re- 
lated to the properties of motor proteinsi^"^. 

Without permanent crosslinks, the filament network can 
plastically evolve and macroscopic flow may emerge, result- 
ing in non-trivial pattern formation as observed in in vitro sys- 
tems consisting of microtubules and various associated mo- 
tors-^"—, and actin-myosin complexes-!*- This richer problem 
has inspired the development of analytical theories incorpo- 
rating filament flow and active driving^"—. One approach 



This journal is ©The Royal Society of Chemistry [year] 



Journal Name, 2010, [vol], 1-Q4] |1 



extends the hydrodynamic equations of liquid crystals in their 
nematic phase- 8 to include active terms (with phenomenologi- 
cal coefficients) obeying the required symmetries^ - — . These 
active nematodynamic equations admit stable solutions with 
line defects in the director field, include cylindrical spirals 
that rotate due to active processes^; these were likened to 
the vortices seen in quasi-2D microtubule experiment s 16 ' 17 , 
although no quantitative comparison has yet been made. An 
alternative approach extends the Smoluchoswki equations for 
rigid rods— to include active driving— - —, and although these 
active terms were derived from microscopic considerations, 
some coarse-graining is still required, resulting in what can be 
regarded as a mesoscale model. 

What is lacking despite this plethora of modelling is a 
clear picture of the microscopic mechanisms driving the self- 
organization observed in experiments. Ideally the large-scale 
equations could be found by coarse-graining a suitable micro- 
scopic model, but this would inevitably involve approxima- 
tions whose validity would need to be assessed. Furthermore, 
the only coarse-graining attempted so far failed to generate all 
of the terms expected on symmetry grounds, with implications 
for the formation of vortices^. It is into this state of affairs 
that simulations can play a key role, permitting as they do full 
control and access of all microscopic details. Non-Brownian 
simulations mimicking the microtubule experiments generate 
asters and apparently vortice s 16 i 17 i 4Q i 41 , but did not include ex- 
cluded volume interactions between filaments and thus lacked 
nematic elasticity, making comparison to the nematodynamic 
theories^"— problematic. Further simplified models with- 
out filament growth found no defects^. Point-like defects 
were also claimed in two-dimensional simulations of inelastic 
rods^"^, but these were somewhat coarse-grained and again 
provided no microscopic picture for the role played by indi- 
vidual motors. 

Here we present the results of extensive simulations of a 
two-dimensional model for active gels, with the goal of com- 
plementing existing analytical and experimental approaches 
by providing a microscopic underpinning for any observed 
non-trivial macroscopic pattern formation. Our chosen numer- 
ical scheme includes: (i) Mobile motors as the originator of all 
non-equilibrium effects; (ii) Excluded volume interactions be- 
tween filaments, so nematic elasticity is present; (iii) Thermal 
diffusion of filaments, which will be relevant to actin if not 
microtubule systems; (iv) Mechanical elasticity of filaments, 
which can store and release elastic energy in the form of bend- 
ing and stretching modes. Filaments are polar, that is they 
have well-defined [+] and [— ]-ends, and motors move strictly 
towards [+]-ends, breaking microscopic reversibility. Some 
snapshots are given in Fig. Q] The structures evident in these 
snapshots, and associated dynamic anomalies, are fully char- 
acterized below. Given the large aspect ratio of single fila- 
ments, finite-size effects are pronounced and we have taken 



great efforts to control variation with system size for all of the 
quantities considered. 

One of our central findings is that the action of the motors 
can drive the formation of non-trivial structure formation in 
one region of parameter space, and anomalous dynamics in a 
different, albeit overlapping parameter regime. Motor-driven 
structure formation is most evident for high motor densities, 
where the nematic ordering breaks down and filament clus- 
ters form. The nature of these clusters depends on the motor 
speed. For motors sufficiently fast to dominate over thermal 
diffusion, localized clusters form in which the filaments are 
bound at their [+]-ends, as evident in Fig. |TJb). This can be 
viewed as non-equilibrium polar ordering on scales compara- 
ble to the filament length L, that becomes isotropic on much 
larger lengths. The dynamics of such states is anomalous in 
that the mean-squared displacement of filament centers ex- 
hibits a superdiffusive regime in which it scales faster than lin- 
early in time, but this effect becomes less pronounced as the 
motor speed increases, a non-intuituve finding that we explain 
below in terms of a vanishing population of mobile motors. 

Maintaining a high motor density but decreasing the motor 
speed so that motor motion and thermal effects compete, we 
find a region of parameter space that exhibits both anomalous 
diffusion and non-trivial structure formation; see Fig. Old). 
Extended, polarised clusters form that can span the system 
size, generating large voids between the high-density clusters 
that is evident as an increase, possibly a divergence, of the 
static structure factor S(q) as the wavenumber q — > 0. Su- 
perdiffusive mass transport is most pronounced in this param- 
eter regime. Nonetheless the polarity of filaments separated by 
distances comparable to a filament length or more are uncor- 
rected, in contrast to the prevalent assumption of continuum 
modelling that the polarity field is slowly-varying over such 
lengths. 

We also identify a third region of parameter space in which 
the system is structurally similar to the equilibrium limit with 
no motors, exhibiting nematic ordering similar to Fig. [TJ e X 
but nonetheless displays anomalous, superdiffusive transport. 
This regime corresponds to the same motor speeds that gener- 
ate the extended clusters described above, but a lower density 
of motors. The key observation underlying this quandary is 
that although the nematic ordering of the filaments appears 
normal, the polarity of nearby filaments are nonetheless cor- 
related over distances corresponding to a number of filament 
diameters. Since the coupling between motors and filaments 
is polar, such correlations can result in a persistent direction of 
motor forces acting between filaments, generating net transla- 
tional motion that can be viewed as a form of emergent self- 
propulsion. This intuitive idea is supported by showing fil- 
ament velocity and polarity are most strongly correlated for 
parameters for which anomalous transport reaches its maxi- 
mum. 
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This paper is arranged as follows. In Sec. |2] we describe 
in detail the numerical model, before turning into Sec. [3] to 
briefly describe the behavior of passive systems in which mo- 
tors do not move. This provides a comparison for the results 
of active systems presented in the succeeding section, where 
we separately focus on static quantities in Sec. 14.11 and then 
the dynamics in Sec. 14.21 Despite performing a systematic 
investigation of parameter space, and including the physical 
mechanisms thought to be necessary for the emergence of the 
structures seen in the microtubule experiment s 16 ' 17 , we never 
observe vortices. Possible causes are discussed in Sec. [5] 

2 Model 

We are interested in determining the universal features of 
motor-driven filament systems, independent of the atomic de- 
tails of either the filaments or the motors (or motor clusters). 
To this end, we describe both filament segments and motors in 
a simplified manner requiring only a few degrees of freedom, 
and focus on the collective effects of many motors and fila- 
ments. See Fig.|2]for a schematic representation of essential 
aspects of the model, which we now describe in detail. 

Filaments are described as polar semi-flexible polymers of 
monomers with lattice spacing b. Rigidity and contour length 
are maintained by elastic energies penalizing local fluctuations 
in both monomer separation and bending. For the former, two 
monomers instantaneously separated by a distance r incurs a 
cost 5£ stretch = %a 2 E(r-b) 2 /2b, where E is the material's 
Young's modulus and a the radius of the cross-section (as- 
sumed circular)^. In addition, the energy for bending by a 
local angle (]) is given by 8£' bend = K$ 2 /2b with K the bend- 
ing modulus, which for a rod with a circular cross-section of 
radius a is given by K = nEa 4 /4^. Once thermal fluctua- 
tions from the solvent are included (see below), and assum- 
ing sufficient resistance to stretching to inhibit contour length 
fluctuations, this gives a persistence length i p « K/lk-gT in 2 
dimensions. For our parameters, i p /L ps 10/3 in terms of the 
filament length L, so our filaments are semiflexible. 

Excluded volume is incorporated as a repulsive potential be- 
tween non-bonded monomers. We use the repulsive part of the 
Lennard-Jones potential with range a and energy e, truncated 
and shifted to ensure continuity of the first derivative at the 
maximumrange 2 1 / 6 o, i.e. V ev (r) = 4e(o/r) 6 [{o/rf - l] + 
e, with r the distance between monomer centers. To avoid a 
proliferation of parameters we simply set a = b, while £ = 
Sk^T is set sufficiently high to ensure excluded volume dom- 
inates over thermal fluctuations for overlapping monomers. 

Thermal fluctuations arise due to the interaction between 
the filaments and the solvent. Assuming the usual low- 
Reynolds number limit for biophysical systems, we adopt a 
Brownian dynamics algorithm that describes forces due to sol- 
vent fluctuations in the non-inertial limit^-. To better model 




Fig. 1 (Color) Snapshots for systems with varying motor density 
and speeds. The left-hand column (a), (c) and (e) corresponds to a 
low motor density, and the right-hand column (b), (d) and (f) to high 
motor density. Similarly, the top row (a), (b) corresponds to fast 
motors, the middle row (c), (d) to slow, processive motors, and the 
bottom row (e), (f) to static motors. Filaments are shaded light 
(dark) towards their [+] ([— ]) ends, respectively, and the blue-green 
hue is arbitrary. Motors are shown as red lines. The system size was 
4L x 4L in all cases. Larger snapshots and movies are available from 
Ref.— . In terms of the model parameters (see Sec.O, £aAd = 20 
for (a), (c), (e) and 40 for (b), (d), (f); k M x b = 75 for (a), (b), 0.75 
for (c), (d) and for (e), (f). 
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the anisotropic friction coefficients for elongated particles , 
we first project forces parallel and perpendicular to the local 
filament axis, and then apply the usual displacement incre- 
ments with solvent friction coefficients yll and y x = 2y", re- 
spectively. Note that solvent forces and drag obey the usual 
fluctuation-dissipation relation via the temperature kgT, so the 
solvent is in equilibrium; all non-equilibrium effects derive 
from the motion of motors, which we now describe. 

Motors (or motor clusters) are modeled as two-headed 
harmonic springs that are simultaneously attached to two 
monomers; motors with one or no attached heads are not 
explicitly represented. The spring extension is defined in 
terms of the separation between the monomer centers relative 
to its natural length, here taken to be the excluded volume 
range 2 1//6 o, and the spring stiffness is fixed at ksT/b 2 . At- 
tachment, detachment and motor motion are defined by the 
following rate-based rules, (i) Motors attach at a rate Ica be- 
tween any two monomers on different filaments whose centers 
are within a specified distance, which for simplicity we take 
to be the spring's natural length 2'/ 6 a. (ii) Both motor heads 
simultaneously detach at a rate kry that does not depend on 
the positions of heads along the filaments, nor on the spring 
tension, except in that severely stretched motors with elastic 
energies exceeding w 7.5kgT detach immediately. Increased 
detachment rates at filament ends, thought to be important for 
the formation of vortices (but not asters)ii, could easily be in- 
cluded at a later stage, (iii) Each motor head moves in a Monte 
Carlo-like manner: The change in motor spring energy AE for 
a head to move m > 1 monomers towards the attached fila- 
ment's [+]-end is calculated, and this trial move is accepted 
with probability kyie^^l^ per unit time if AE > 0, or £m 
if AE < 0. Each m is drawn uniformly from the fixed range 
[l, OT max], and the acceptance probability is suitably normal- 
ized to ensure invariance with respect to m max . Motors already 
in tension are less likely to move due to the increase in strain 
energy, giving an effective stall force of the order of k%T jb. 
Motors already at [+]-ends simply do not move, but may de- 
tach at the usual rate. Motors do not move if to do so would 
exceed the maximum spring energy described in (ii). 

All filaments have M = 30 monomers and hence are of 
length L = 30b, giving an aspect ratio L/b = 30. Typically 
we use a 4L x 4L box, but to check for convergence with sys- 
tem size for any given quantity, we also simulated 6L x 6L 
and 8L x 8L boxes for a representative sample of parameter 
space, corresponding to 2 vertical and 2 horizontal lines in the 
parameter space discussed below. The area fraction of fila- 
ments was fixed at (j) w 21% throughout, which was checked 
to give a nematic order parameter close to unity in the absence 
of motors. At t = filaments are placed in a smectic con- 
figuration to ensure no significant initial overlap, with each 
filament's polarity independently chosen to be ±p with equal 
probability. Steady state is identified as when the two-time 




Fig. 2 (Color online) Schematic of model parameters. Filament 
polarity is defined by [+] and [— ]-ends as denoted in the figure. 
Motors are 2-headed springs defined by attach, detach and motion 
rates k£> and kw, respectively, where the movement rate is 
attenuated by an exponential factor depending on the increase in 
spring energy AE for the proposed move (AE = if this change is 
negative). See text for details. 

mean squared displacement of filament centres-of-mass x(f ), 
(Ar 2 (f w ,f w +t)} = (|x(f w + t) -x(f w )| 2 ), ceases to depend on 
f w and only varies with the lag time t, i.e. becomes {Ar 2 (t )). 

The results from the simulations will be described in terms 
of the two dimensionless parameters &a/&d and kyftb, where 
Zi, = Lby/Ak^T is the time for an isolated filament's centre- 
of-mass to translationally diffuse over one monomer diam- 
eter. Thus k^/kry gives a rough measure of motor den- 
sity, and kyi%b measures the competing effects of motor mo- 
tion to thermal diffusion over lengths of the monomer spac- 
ing b. We fix the value of k\) to be small relative to xT , 
i.e. koZ/, = 7.5 x 10~ 3 <C 1, and consider £m varying from 
<C kry to kyi ^> kry, where the former case corresponds to 
non-processive motors that detach before significant motion, 
and the latter to processive motors which may move many 
monomers before detaching. 

3 Passive systems: &m = 

We first summarize the behavior of passive systems with 
= 0, before turning to active systems with kyi > 0. With- 
out motor motion, filament polarity is irrelevant and the time- 
averaged effect of motor attachment and detachment is to 
generate an effective attraction between filaments. As the 
ratio &a/£d between attachment and detachment rates in- 
creases, the magnitude of motor-mediated attraction increases 
and induce a heterogeneous density distribution as evident in 
Fig- E 1 )- Increasing k^lk-Q even further gives a percolating 
network that nonetheless does not fully phase separate over 
our simulation time. 

To identify the length scales associated with the den- 
sity fluctuations, we calculate the static structure factor 
S(q) of filament centers, defined as the angular average 
of 5(q) =A^ _1 (|«(q)| 2 ) over q = q/q with q = |q|. Here, 
n(q) = Yla=\ e ' q x is the Fourier-transformed distribution of 
the filament center-of-masses x a , a = 1 . . .N. Results for 
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Fig. 3 (Color online) Static structure factor S(q) for = and the 
&aAd given in the key on log-log axes. The horizontal and vertical 
dotted lines correspond to S(q) = 1 and q = 2k/ L, respectively, with 
L the filament length. Error bars give the spread between 
independent runs. In all cases the system size was 6L x 6L. 
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Fig. 4 (Color online) Mean squared displacement (Ar 2 ) for kyi = 
and A'a/A'd given in the key for the same systems as Fig. [5] The thick 
dashed line has slope 1 corresponding to normal diffusive scaling 
(Ar 2 ) <x t, and the solid horizontal line corresponds to displacements 
equal to the filament length. 



1 < £a/&d < 80 are shown in Fig. [3] For low AaAd, S(q) 
exhibits a small peak around the wavenumber corresponding 
to the filament length L, suggesting a slight degree of smectic 
ordering with filaments aligned end-to-end, before decaying 
quadratically for smaller q. Increasing &aAd to around 40, 
corresponding to o(10) motors per filament, dramatically en- 
hances the height of this peak. Snapshots such as Fig. |TJf) 
reveal local filament bundles that align end-to-end, explaining 
this peak. For £aAd in the range 60 — 80, the S(q) collapse 
onto a single curve with S(q) ~ \ — 3 as (7— >0 and snapshots 
reveal similar pictures of percolating networks. It is clear that 
the system does not phase separate for k^/ko > 60 on our 
simulation times, but instead undergoes kinetic arrest into a 
long-lived metastable gel. We therefore restrict attention to 
^a Ad < 40 for the active systems in Sec. |4] 

One of our key findings for active systems is the existence 
of enhanced diffusion, so for comparison we present in Fig. [4] 
the mean-squared displacements for passive systems. Given 
the lack of motor-mediated driving, the sole effect of motors 
is to bind filaments and thus reduce self-diffusion, and this 
trend is immediately apparent from the figure. For low mo- 
tor densities the mass transport is roughly diffusive, with a 
diffusion constant that decreases with increasing £a/^d and 
the binding between filaments increases. At high attach rates 
the mass transport is substantially reduced, and never becomes 
fully diffusive over our data window, instead giving weak sub- 
diffusion (Ar 2 ) ~ t a with a « 0.9 at the maximum time lags 
available. We nonetheless expect normal diffusion with a = 1 
to be recovered at late times. 



4 Active systems: £m > 

The movement of motor heads along the polar filaments gen- 
erates equal-and-opposite forces that drive relative motion be- 
tween filaments or filament clusters. Non-trivial flows and 
pattern formation may thus become stable. In this section 
we describe the structure and dynamics of active systems in 
terms of the 2 dimensionless parameters &a Ad, which broadly 
corresponds to the density of motors, and the bare motor 
speed kuib- 

4.1 Statics 

Snapshots for systems with mobile motors reveal similar fila- 
ment ordering as for static motors when &aAd remains below 
some crossover value. This crossover value depends on pa- 
rameters such as filament length and motor spring stiffness, 
and for the systems considered here occurs around k^/ku « 
15. When &aAd exceeds this value, the increased motor den- 
sity induces filaments clustering and non-trivial structure for- 
mation, the nature of which depends on the speed of the mo- 
tors. For £m <C ^d, motors detach at a faster rate than moving 
and become non-processive. Filaments form apolar bundles 
similar to the passive case ku = 0. Conversely, for kyiib 3> 1 
motor motion dominates over thermal diffusion and we ob- 
serve filament clusters bound at their [+]-ends as in Fig. [Tib). 
For the intermediate regime k^Zb <C ^mX^ <C 1, motors are 
processive but compete with diffusion in structure formation. 
This results in the extended clusters evident in Fig.[TJd). 

Underlying the observed clustering is a monotonic increase 
of the density of attached motors as the ratio &aAd is in- 
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Fig. 5 (Color online) Number of motors per filament ny^/N versus 
Ad f°r the dimensionless motor rates kyili, given in the key. The 
thick black line segment has slope 1. For all lines, symbols are 
larger than the error bars. 

creased. As plotted in Fig. |5j the number of motors per 
filament is roughly proportional to k^/ko, with a prefactor 
that decreases with motor speed ky\_- This linear dependency 
on £a/&d can be understood as the steady state solution of 
the simple rate equation d t n mot <^ k& — « mo t^D governing the 
number n mot of attached motors between two monomers held 
within the proscribed attachment range. Deviations from strict 
linearity are evident at high £a/^d > o(10) when motors in- 
duce clustering and the global attachment rate becomes cou- 
pled to structure formation. The decrease in motor density for 
high kM is likely due to the greater fraction of motors stretched 
close to the detachment limit mentioned in Sec. [2] but this ef- 
fect is reduced for higher &aAd when the densities become 
similar for all motor speeds. 

4.1.1 Cluster size distributions: To quantify the cluster- 
ing apparent in snapshots, we first consider the distribution 
P(n c ) of clusters consisting of n c filaments, averaged over 
time and independent runs, where two filaments are consid- 
ered to belong to the same cluster if there is at least one mo- 
tor simultaneously attached to both. The mean cluster size 
(»c) = Hn c n cP{nc) is plotted in Fig.[6j and reveals a mono- 
tonic increase with kp^jk-a following a roughly exponential re- 
lationship. It also decreases with increasing motor speed ky\, 
in part because the density of motors decreases (see Fig. [5}, 
and also because they drive the formation of localised fila- 
ment clusters bound at their ends. Increasing the system size 
increases cluster sizes for high &aAd, but the effect is small 
on logarithmic scales. 

Considering only the mean cluster size (n c ) obscures the 
fact that the shape of the full distribution P(n c ) qualitatively 




Fig. 6 (Color online) Mean number of filaments in a cluster (n c ) 
versus &aAd f° r the motor speeds given in the key. The system size 
was fixed at 4L x 4L and error bars are smaller than the symbols. 

changes as the parameter space is traversed, see Fig. UJ For 
most of the parameter space considered, P(n c ) is well approx- 
imated by a simple exponential decay, and does not vary with 
system size. However, for high £aAd an d k^b ~ 1> the dis- 
tribution P(n c ) deviates from an exponential, in two ways. 
Firstly, a second peak corresponding to system-spanning clus- 
ters with « c « N emerges for ku^b ~ 1 ■ Secondly, and only 
for slow, processive motors with koib ~ ku%b ~ 1, the small-n c 
decay becomes power law rather than exponential. Although 
noisy, this algebraic decay can be approximately fitted by an 
exponent -2, as shown in Fig. [7] We note that the region of pa- 
rameter space for which P(n c ) exhibits power-law decay co- 
incides with those parameters that exhibit anomalous small- 
wavenumber density fluctuations, as described in Sec. 14.1.21 

The magnitude of the second peak decays monotonically 
with increasing system size, suggesting it will vanish in the 
limit of large system size. This is apparent in plots of the inte- 
grated area of the second peak, P(n c >N/2) — L^ =JV /2^ > ( Hc )' 
which can be understood as the probability to find a cluster 
that is comparable in extent to the system size. The lower 
cut-off N/2 is arbitrary, but as long as it falls into the middle 
region where P(n c ) ~ 0, as here, its precise choice does not 
measurably alter the result. As shown in the inset to Fig. |7j 
this quantity decays roughly as A/ -1 . P(n c > N/2) also mono- 
tonically decreases with increasing kyflb as show in the figure, 
before vanishing entirely when ku^b > 0(1). 

4.1.2 Small-wave number density fluctuations: Al- 
though useful for considering connectivity, the cluster distri- 
bution P(n c ) gives no information regarding spatial distribu- 
tion of clusters, and in particular exhibits no signature corre- 
sponding to the large void formation evident in snapshot of 
Fig- Hd) for high £aAd an d k^ij, ~ kyi%b ~ 1- To consider 
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Fig. 7 (Color online) Log-log plot of the distribution P(n c ) of 
cluster sizes n c for &aAd = 40, k^Z/, = 0.075 and the system sizes 
given in the key. The thick, dashed diagonal line has a slope of —2, 
and the vertical lines correspond to n c = N with N = 175, 400 and 
711 the total number of filaments in the system. (Inset) P(n c > N/2) 
versus N on log-log axes for (from top to bottom) ky^b = 0, 
7.5 x 10~ 3 , 7.5 x 10~ 2 and 7.5 x 10" 1 . P{n c > N/2) = for higher 
kyfii, . The thick black line segment has a slope of — 1 . 



the distribution of mass centers it is useful to look at the static 
structure factor S(q) as denned in Sec. [3] Fig.[8]shows the vari- 
ation of S(q) with kyi^b for fixed &a/£d = 40. Focussing on 
the S(q — > 0) behavior reveals a non-monotonic dependence 
on motor speed: for fast motors kyfib > 0(1), or slow, non- 
processive motors k^ < fcpj, S(q) decays to some small value 
as q — > 0. For intermediate ku, S(q) increases, possibly di- 
vergently, with decreasing q, which coincides with the void 
formation evident in snapshots. 

The limited number of data points and poor statistics for 
small q makes it difficult to characterize the behavior of S(q) 
in this limit. Nonetheless some semi-quantitative observa- 
tions can be made by fitting each S(q) to the form A + Bq c + 
(j e -{g/[2n/L]-l) /D oyer tne ran g e a < 2q-^, which provides a 
stable fit for all parameters considered. Although the fitted 
exponent c is somewhat susceptible to finite size effects, we 
consistently observe values close to c = 2 for low motor den- 
sities £a/^d ~ 1, or very slow or fast motors, £m -C £d and 
kyitb ^ 1 respectively. For high k^/k-Q and intermediate ^m^, 
c decreases and becomes negative. Although the magnitude of 
the fitted exponent depends strongly on system size, the sign 
does not, suggesting some form of small-wavenumber struc- 
ture will persist for large system size. 

A divergent S(q) was predicted from hydrodynamic equa- 
tions of active nematic systems with an exponent c = —2^, 
and was interpreted as emergent directed motion deriving 
from the coupling of elastic modes with the breaking of time- 
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Fig. 8 (Color online) Static structure factor S(q) for Aa/&d = 40 
and the motor speeds given in the key, for system sizes 8L x 8L. For 
clarity only error bars for k^ti, = 0.075 are shown; others are 
comparable. The horizontal and vertical dotted lines corresponds to 
S(q) = 1 and q = 2%/L respectively. 

reversal invariance in driven states^. We shall later argue for 
a form of emergent directed motion in our system when dis- 
cussing the dynamics in Sec. 14.21 which derives however from 
the polar coupling between motors and filaments. It is possi- 
ble that emergent directed motion drives the small wavenum- 
ber density fluctuations in both systems. As described in 
Sec. 14.2.11 the parameter values for which exponents c < 
arise coincide with those for which superdiffusive mass trans- 
port is most pronounced. It is interesting to note that a causal 
link between superdiffusion and small wavenumber fluctua- 
tions has been hypothesized from systems lacking orienta- 
tional degrees of freedom—, in broad agreement with our find- 
ings. 

4.1.3 Nematic and polar ordering: We now move be- 
yond density fluctuations and consider the orientational de- 
grees of freedom. Let p" denote the unit polarity vector from 
the [— ]-end to the [+]-end for each of the N filaments a = 
I ...N. From this can be defined the traceless nematic order 
tensor Qfj = pfpj — The two-dimensional nematic order 

parameter S2D 1S tnen defined by (S2d) 2 = 2 II (Qij) || 2 , where 
(Qij) = N^^aQfj an d ||---l| 2 denotes the Euclidean norm 
||A,j || = £,yAy. S2D varies from to 1, with 1 for perfect 
nematic ordering and when there is no net preferred orienta- 
tion, such as for an isotropic state, asters etc. A contour plot of 
S2D for kp^/kv and ku% is given in Fig. [9] and confirms what 
is visible from snapshots in Fig.Q] namely that nematic order- 
ing breaks down for &a/£d ~ 15 and processive motors with 
kyi 3> kx). This data is for the smallest system size, for which 
we have most data points; the effects of system size have been 
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checked and confirms the crossover to low S2D is robust. The 
value of S2D may be approaching zero for points in the upper 
right of this plot, but the poor statistics and limited number of 
points rules out extrapolation to the infinite system-size limit. 

The nematic order parameter 620 tell us nothing about cor- 
relations in the polarity of filaments. To consider correla- 
tions in filament polarity, we must go beyond this nematic de- 
scription and consider order parameters that respect reverses 
in polarity p — > —p. Firstly we note that the mean polarity 
(p) = iV _1 LaP a vanishes for the entire parameter space con- 
sidered here. This is in agreement with the theoretical predic- 
tion that states with non-zero mean filament polarity are only 
stable if the motors are polar, i.e. have attachment rates that 
depends on the relative filament polarity--, unlike the apolar 
motors considered here for which the attachment rate depends 
only on separation. 

It is natural to consider spatial correlations in polarity along 
directions parallel to the filament axis separately to those in 
perpendicular directions. We therefore define the perpendicu- 
lar polarity correlation function Cp p {r) as 

c± (r) Ia,pp"-p p 5(Mx a ~xP|)sin 2 9 
PP L a , p S(r-|x«-xP|)sin 2 e 

(1) 

where is the angle between p 01 and xf 1 — x 01 (pP could equally 
have been chosen due to the symmetry of Eq. ([TJ). The corre- 
sponding parallel function C pp (r) is defined analogously, with 
the sin 2 weighting factors replaced by cos 2 8. Both pro- 
jections exhibit approximate exponential decay with r, with 
longer-range correlations for higher motor densities &a/&d, as 
demonstrated in Fig. [10] These correlations become shorter 
range as kyi — > 0, confirming the smooth approach to the pas- 
sive limit kyi = when all polarity correlations trivially van- 
ish. It should be noted that significant polarity correlations are 
not mutually exclusive with nematic ordering. With regards 
to finite-size effects, the data shows no apparent trend with 
system size for all of the parameter space except for the high 
£a/£d and k^Tb ~ k^Xj, ~ 1 when extended polar clusters arise, 
for which the range of polarity correlations was still increas- 
ing for the largest system sizes simulated. Finally, we never 
observe significant correlations on lengths larger than the fila- 
ment length L, which will be discussed in Sec.|5]in relation to 
hydrodynamic modelling. 

4.2 Dynamics 

Motor motion represents an energy flux driving relative fila- 
ment motion in violation of the principles of thermodynamic 
equilibrium, and might therefore be expected to result in 
anomalous transport properties. Data confirming this expec- 
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Fig. 9 (Color online) Nematic ordering parameter Sjd f° r ^aAd 
and ky[Zj, for fixed system size 4L x 4L. In this figure, and in Fig. 1121 
below, the white discs denote actual data points; linear interpolation 
within the enclosing triangle is used between points. The thick 
horizontal dashed white line corresponds to ky[ = &d- 




Fig. 10 (Color online) The transverse polarity correlation function 
Cp p (r) for kyftb = 0.075, system size 4L x 4L and the attachment 
rates given in the key. For clarity only errors bars for &aAd = 40 
are shown; the others are comparable or smaller. The thick dashed 
line is proportional to e^''/^/ 2 ! . 
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Fig. 11 (Color online) Emergence of directional motion on the 
one-filament level under the action of motors, denoted by the short 
horizontal lines, (a) Net zero displacement for parallel filaments, 
(b) Net displacement towards the [— ]-end for both filaments in an 
anti-parallel pairing. 



tation is provided in this section. The key observation under- 
lying these findings is that the coupling between motors and 
filaments is polar, in that motor heads move strictly towards 
filament [+]-ends. Filaments are thus expected to move in the 
direction of their [— ]-ends in response to the motor forces, as 
seen by considering connected parallel and anti-parallel fila- 
ment pairs as in Fig. [TT] This effect also emerges from one- 
dimensional continuum modeling 54 ' 55 . Thus filaments may 
move persistently in a direction that is coupled with their po- 
larity. This scenario is reminiscent of self-propelled parti- 
cles that are known to typically exhibit enhanced mass diffu- 
sion, such as an anomalous scaling with time of displacements 
transverse to particle motion^"— or a longitudinal diffusion 
constant increased by the activity—. In this light, anomalous 
transport should also be expected here. 

To confirm the correlations between filament motion and 
polarity, the filament velocity must be defined, for which it 
is necessary to average over a finite time interval since there 
is no instantaneous velocity in Brownian dynamics. We de- 
fine the velocity of filament a at time t to be v" = [x a (t + 
At) — x a (f)]/Af where At m 65%b, corresponding to trajecto- 
ries spanning a few monomer diameters. Correlations between 
polarity and velocity are then easily quantified as (v • p) /v 
with v = \/(v ■ v) the mean filament speed. This is plotted 
in Fig. [T2] and shows higher correlations for large kx/ku 
and kyi%b ~ 10 — 10 _1 , which we will show below is also 
roughly where the degree of anomalous mass transport is high- 
est. Note also that the correlations are negative, confirming 
filaments move in the direction of their [— ]-ends. 



^'aA'd 

Fig. 12 (Color online) Correlation between individual filament 
velocity and its polarity, (v • p) /v, versus A;a/^d an d ^M'tfe- As this 
quantity is relatively insensitive to system size, the largest available 
system was selected for each point. 



4.2.1 Super-diffusive mass transport: Mass transport is 
conveniently quantified by the mean squared displacements 
of filament centers, (Ar 2 (t)}, where for simplicity we average 
over filament polarity. For the available time window, which 
typically corresponds to displacements from a fraction of b 
to a few L, a slow, diffusive or sub-diffusive region is ob- 
served at small times, sometimes followed by a crossover to 
super-diffusive scaling at later times in which (Ar 2 (f)) ~ t a 
with a > 1. See Fig.[l3]for some examples. The variation of 
these curves with system size is weak or non-existent. As with 
the passive case in Sec. [3] we presume an eventual crossover 
to normal diffusion with a = 1 for lag times exceeding our 
achievable time window. 

To quantify the extent of deviation from normal diffusion 
as a function of the microscopic parameters, it is convenient 
to condense the variation of logarithmic slope a over the data 
window into a single scalar. To do this, we first smooth the 
data by fitting each curve to the sum of two power laws, 
(Ar 2 ) = Cit c ' +C2t C2 , which gives a reasonable fit in all 
cases. We then take the logarithmic slope of this fit at the 
point when (Ar 2 ) = L 2 . The result is given in Fig. [14] and 
shows a decrease in anomalous diffusion away from the peak 
value at high k^/ko and motor speeds roughly in the range 
kuXb ~ kyfib ~ 1 ■ Choosing a different point along the curve to 
extract the slope results in different values but similar trends. 

A striking feature of Fig.[l4]is the non-monotonicity of the 
degree of anomalous diffusion with motor speed kyi, which 
approaches normal diffusion for ky[if, 3> 1 ■ The reason is not 
hard to find. As evident in Fig.[TJb), when motor motion dom- 
inates over thermal diffusion, they rapidly reach the filament's 
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Fig. 13 (Color online) Mean squared displacements divided by 
time, (Ar 2 )/t, so that normal diffusive scaling (Ar 2 ) oc / shows up as 
a horizontal line. Here k^Zf, = 0.075, the system size was 8L x 8L 
and the Ad are given in the legend. The thick dashed lines are 
fits to C\t Cl +C2t C2 for each set of data points. The black horizontal 
line is to guide the eye, whereas the diagonal black line corresponds 
(on these axes) to (Ar 2 ) = L 2 and confirms trajectories exceed the 
filament length for these data sets. 




Fig. 14 (Color online) Scaling of the mean squared displacement 
with time t, i.e. a = dln(Ar 2 (r)) /3lnf, at the point when 
(Ar 2 ) = L 2 , versus f° r me Ad given in the key. For clarity 
only error bars for #aAd = 20 are shown; others are comparable 



[+]-end where they then stall. Indeed, the fraction of motors 
with at least one head attached to a [+]-end never drops be- 
low 80% for kyi%b > 1. Since such heads can no longer move, 
they do not drive filament separation and the net motor activity 
decreases, restoring normal diffusion. 



4.2.2 Velocity correlations: Using the same definition of 
filament velocity v a as above, it is possible to calculate spatial 
correlations in velocities projected parallel c|,,(r) and perpen- 
dicular C^,(r) to the filament axis. Here, c|A(r) are defined 
analogously to the polarity correlations in Eq. (Q]), except with 
the p" replaced by v a /v with v the mean filament speed as 
before. Examples are given in Fig. [15] demonstrating a grow- 
ing range of correlated motion as the motor density increases. 
This trend was observed throughout the parameter space con- 
sidered. 

Also shown Fig.[T3]is an example of the variation with sys- 
tem size for the highest £a/&d, which demonstrates conver- 
gence for the two largest system sizes considered. The decay 
is exponential, approximately of the form sr r l\ L l^ and thus 
becomes negligibly small on lengths on the order of the fil- 
ament length L. We never observe correlations decaying on 
lengths larger than L, and conclude velocities on lengths much 
larger than the filament length are uncorrelated, a point that is 
discussed in Sec. [5] 
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Fig. 15 (Color online) Velocity correlations C^,(r) perpendicular to 
the filament axis, for k^Z/, = 0.075 and the £aAd ar, d system sizes 
given in the legend. The thick dashed line is °c e~ r// [ L / 6 ]. 
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Fig. 16 (Color online) Example of asters at low density with area 
fraction <|> ~ 5.3%. The system parameters were A'a/^d = 250, 
ktA^b = 75. Motors are short straight lines concentrated at filament 
[+]-ends. 

5 Discussion 

It is apparent from the results described above that for the con- 
sidered density of (]) ps 21%, structures similar to the asters 
and vortices in the microtubule experiment s 16 i 17 are never re- 
produced. The absence of asters may be a simple matter of 
density. The end-bound clusters formed by rapid motors in 
Fig- Eh) are not able to form circular structures due to steric 
hinderance with nearby clusters. Lowering the filament den- 
sity removes this effect, permitting full asters consisting of a 
single layer of filaments to form, as demonstrated in Fig. [16] 
However, vortices did not arise at lower densities, even when 
increased motor detachment rates at [+]-ends was included. 
We note that a continuum model— that extended an earlier 
version 3 - to include, amongst other features, a form of steric 
hinderance between filaments, favored asters over vortices rel- 
ative to the earlier work, suggesting excluded volume may also 
be a factor. Here, we expand upon the relationship between 
our numerical findings and the related theory and experiments 
that has already been touched upon, with an emphasis on pos- 
sible reasons for the non-observation of vortices. 

5.1 Comparison to experiments 

In their in vitro actin-myosin experiments, Backouche et al. 
claimed that active structure formation was only possible with 
the addition of a small concentration of the passive crosslinker 



fascin—. A small fraction of passive crosslinks between adja- 
cent filaments has also been predicted to significantly increase 
their rate of alignment by mobile motors 4 **. Since no passive 
crosslinkers were included in our simulation, this might be a 
factor in the qualitatively different structures found here. Fur- 
thermore, fascin is a bundling protein that promotes the forma- 
tion of polar actin bundle s 64 ! 65 , thus its dominant role may in 
fact be to permit non-zero mean polarity. If so, it may play the 
role of the polar motor clusters that preferentially bind to par- 
allel (as opposed to anti-parallel) filaments, which were shown 
theoretically to stabilize homogenous polarized states and pro- 
duce a correspondingly richer stability diagram—. However, 
there is no reason to suppose that polarity-dependent binding 
of active or passive components played a role in the micro- 
tubule experiments, and was not included in the associated 
simulations 16 ^ 7 -, suggesting this is not a reason for the non- 
emergence of vortices in our work. 

Apart from molecular motors, a second form of non- 
equilibrium activity in biofilament gels is the spontaneous 
translation of filament center-of-mass via unequal monomer 
addition rates at opposite ends, known as treadmilling and 
prevalent in vivc&. If present, this would place the system into 
the class of self-propelled particles, for which a richer variety 
of structure and dynamics is expected^ 5 - 6 - " 59 ' 66 ' 67 . However, 
treadmilling was thought not to play a role in the actin-myosin 
experiments™, and was inhibited to some extent by taxol in 
the microtubule experiments (although microtubule growth 
was present) 16 ' 17 . Until the role of treadmilling or filament 
growth is categorically denied experimentally this remains a 
possible missing factor, but without further experimental guid- 
ance we can merely speculate on its possible role. Including 
such features numerically should be straightforward, and in- 
deed has already been performed for simulations of branching 
networks^"—. 

Finally, the role of dimensionality should not be over- 
looked. Our simulations are strictly two-dimensional, whereas 
the actin-mysoin experiments were three-dimensional— and 
the microtubule experiments were performed in a thin 5/jm 
chamber that enforces almost-parallel filament orientation 
while facilitating overlap, which can be regarded as a 
quasi-2D syste m 16,17 (the associated simulations were two- 
dimensional but without excluded volume, thus also permit- 
ting free overlap). The interaction terms in the hydrodynamic 
and mesoscale models were based on a three-dimensional ker- 
ne l 38 ' 39 . Thus the excluded volume interactions included in 
our simulations may be far more strict than the experiments or 
models, possibly having an inhibitory effect on structure for- 
mation. To numerically probe higher dimensionality is com- 
putationally expensive but should be possible with a restricted 
sampling of parameter space. 
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5.2 Comparison to hydrodynamic models 

An assumption common to many of the analytical models is 
that the velocity and polarity fields are hydrodynamic, in the 
sense that they have long-wavelength components extending 
across lengths much larger than the filament length L. In con- 
trast, as described in Sees. 14. l l and 14.21 above, we never ob- 
serve polarity or velocity correlations that decay on lengths 
larger than L. It is possible that long-wavelength correlations 
emerge at much higher filament densities, but we were not 
able to check this so far due to computational limitations. Al- 
ternatively, our numerical scheme may be oversimplified, in 
that momentum is not conserved by the solvent-filament in- 
teractions. Since the standard argument for the emergence of 
a hydrodynamic velocity field requires momentum conserva- 
tion—, this may explain its absence, but correcting for this 
omission numerically is difficult and will require ad hoc mod- 
ification of the driving terms such as in Ref.— , or the incorpo- 
ration of full hydrodynamic interactions—. 

A slowly varying polarity field is assumed in all theoretical 
models and coarse-graining schemes- 2 ^ - — . The short-ranged 
decay of the polarity correlations in our simulations therefore 
makes the comparison of the results difficult. This may be 
partly a question of time scale. For the density considered 
here, the rotational diffusion time for filaments to significantly 
change orientation can be large. This is particularly true when 
the nematic order parameter S2D is close to unity, when the 
orientation autocorrelation function decays by as little as 5% 
over a production run. This may inhibit the formation of large, 
polarity-correlated regions. Nonetheless even in this region of 
parameter space, motor-driven translational separation of fil- 
aments according to their polarity does occur. Furthermore 
when 5id "C 1, filament rotation is substantially enhanced. 
The lack of long-wavelength polarity correlations is therefore 
somewhat of a mystery. 

The nematodynamic theories take as an input parameter 
the active stress generated on filaments, whose sign deter- 
mines if this stress is contractile or extensil e 2Q i 23 . Our micro- 
scopic approach does not impose the sign of the force dipoles 
generated by motors, so instead it must be measured. The 
mean dipole moment acting between filament pairs is defined 
as K = jyJ— L(ap) F a P ■ (x^ — x a ), where the sum is over all 

Ni a ^ filaments pairs (a, p) connected by at least one mo- 
tor, and F a P is the total force on a due to motors connect- 
ing a and p. Thus K > corresponds to contractile dipoles, 
and K < to extensile ones. Preliminary results are given in 
Fig- El demonstrating uniformly contractile stress K > with 
a magnitude that reaches a maximum for intermediate motor 
speeds k^lb < ku^b < 1- An alternative measure in which 
both F"P and xP - x" are first projected along the filaments' 
axes before summing, as defined in the figure caption, fol- 
lows the same trend. Thus our motor rules lead to contractile 




Fig. 17 (Color online) Variation of the dipole moment with for 
&a Ad = 40 and system size 4L x 4L in steady state. The upper line 
gives the mean dipole moment K = j^L- £(a|3) F a P ■ (xP — x a ) as 
defined in the text, and the lower line gives the projected equivalent 
K proj = ^ £(„p) { (F«P ■ p«)(Ax«P ■ p«) + (FP« ■ pP)(AxN ■ pP)} 
with p a , pP the polarity unit vectors for filaments a and p resp. 

active stresses, theoretically predicted to give the richest be- 
haviour^^. There is also a slight but measurable reduction 
of « 0.4% in the mean filament contour length when K is at its 
highest, confirming contractility. 
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